﻿#define FaceMinCor
//#undef FaceMinCor
using System;
using System.Windows.Media.Media3D;

namespace LxCFD
{
    public struct FvFace
    {
        public FvFace(int[] nodes, int c0, int c1, bool isRightHand)
        {
            NodeIDs = new int[nodes.Length];
            Array.Copy(nodes, NodeIDs, nodes.Length);
            this.OSideCell = c0;
            this.NSideCell = c1;

            this.Center = new Vector3d();
            this.Normal = new Vector3d();
            this.Area = 0;
            this.OToN = new Vector3d();
            this.OToNLength = 0;
            this.gDiff = 0;
#if FaceMinCor
#else
                this._gT = 0;
#endif
            this.gOwner = 0;
            if (!isRightHand)
            {
                Array.Reverse(this.NodeIDs);
            }
        }
        //拓扑信息
        public int[] NodeIDs;
        public int OSideCell;
        public int NSideCell;
        //几何信息
        public Vector3d Center;
        public Vector3d Normal;//lenth为面的面积
        public double Area;
        public Vector3d OToN;//由cell0 中心指向cell1中心的矢量
        public double OToNLength;
        //用于计算
        public double gDiff;
#if FaceMinCor
        public double gT => gDiff;//计算非正交修正时用到
#else
            public double _gT;
            public double gT => _gT;//计算非正交修正时用到
#endif
        public double gOwner;//owner cell 到面的距离 投影到ownerToNeighbor，然后除以ownerToNeighbor的长度
        public void CalGeometryInFace(Point3DCollection points)
        {
            Center = new Vector3d();
            Normal = new Vector3d();
            double areaWeightSum = 0;//如果这几个三角面不在一个平面，不至于出错
            for (int i1 = 2; i1 < NodeIDs.Length; i1++)
            {
                CalTriangleGeometry(out Vector3d CenterMultiply3, out Vector3d NormalMultiply2, points[NodeIDs[0]], points[NodeIDs[i1 - 1]], points[NodeIDs[i1]]);
                double areaWeight = NormalMultiply2.CalLength();
                areaWeightSum += areaWeight;
                Center.AddToSelf(areaWeight * CenterMultiply3);
                Normal.AddToSelf(NormalMultiply2);
            }
            Normal.MultiplyToSelf(0.5);
            Area = Normal.CalLength();
            Center.MultiplyToSelf(1 / (3 * areaWeightSum));
            if (Area == 0)
            {
                LxConsole.Error("存在面积为0的面");
            }
        }
        private void CalTriangleGeometry(out Vector3d CenterMultiply3, out Vector3d AreaMultiply2, Point3D p0, Point3D p1, Point3D p2)
        {
            Vector3d v01 = Vector3d.From2Point3D(p0, p1);
            Vector3d v12 = Vector3d.From2Point3D(p1, p2);
            AreaMultiply2 = Vector3d.CrossProduct(v01, v12);
            CenterMultiply3 = new Vector3d(p0.X + p1.X + p2.X, p0.Y + p1.Y + p2.Y, p0.Z + p1.Z + p2.Z);
        }
        private void calO2N(FvCell[] cells)
        {
            OToN = cells[NSideCell].Center - cells[OSideCell].Center;//nCellCenter - oCellCenter;//
            OToNLength = OToN.CalLength();
            double normalDotOToN = Vector3d.DotProduct(Normal, OToN);
            if (normalDotOToN < 0)
            {
                LxConsole.Error("面的Owner关系与面的朝向不一致，极可能是网格导入错误");
            }
            if (OToNLength == 0)
            {
                LxConsole.Error("网格中心重叠，计算会出现错误！");
            }
            //---------------------修改点1------------------------------------
#if FaceMinCor
            //!!!!!!!!!!!以下为minCorrection
            this.gDiff = normalDotOToN / Math.Pow(this.OToNLength, 2);
#else
            //!!!!!!!!!!!以下为overCorrection修改点1
            this.gDiff = Math.Pow(Area, 2) / normalDotOToN;
            this._gT = normalDotOToN / Math.Pow(this.OToNLength, 2);
#endif
        }
        public void CalGeometryAcrossFace_inner(FvCell[] cells)
        {
            calO2N(cells);
            if (OToNLength == 0) { gOwner = 0.5; return; }//上面已提示过了
            Vector3d ownerToFace = Center - cells[OSideCell].Center;
            gOwner = Vector3d.DotProduct(ownerToFace, OToN) / OToN.CalLengthSquared();
            if (gOwner > 1)
            {
                gOwner = 1;
            }
            else if (gOwner < 0)
            {
                gOwner = 0;
            }
        }
        public void CalGeometryAcrossFace_bound(FvCell[] cells)
        {
            calO2N(cells);
            gOwner = 1;
        }
        public double CalAverage(double[] phi)
        {
            return phi[OSideCell] * (1 - gOwner) + phi[NSideCell] * gOwner;
        }
        public Vector3d CalAverage(Vector3d[] gradPhi)
        {
            return (1 - gOwner) * gradPhi[OSideCell] + gOwner * gradPhi[NSideCell];
        }
        private Vector3d calCFCorrectPart(ref Vector3d aveGrad, double[] phi)//求解C F 网格值修正部分
        {
            return (phi[NSideCell] - phi[OSideCell] - Vector3d.DotProduct(aveGrad, OToN)) / Math.Pow(OToNLength, 2) * OToN;
        }
        public Vector3d CalCFCorrectedGradient_inner(Vector3d[] gradPhi, double[] phi)
        {
            Vector3d aveGrad = CalAverage(gradPhi);
            Vector3d cor = calCFCorrectPart(ref aveGrad, phi);
            return aveGrad + cor;
        }
        public Vector3d CalCFCorrectedGradient_bound(Vector3d[] gradPhi, double[] phi)
        {
            Vector3d aveGrad = gradPhi[OSideCell];
            Vector3d cor = calCFCorrectPart(ref aveGrad, phi);
            return aveGrad + cor;
        }
        public Vector3d CalCFGradCorrection_inner(Vector3d[] gradPhi, double[] phi)
        {
            Vector3d aveGrad = CalAverage(gradPhi);
            return calCFCorrectPart(ref aveGrad, phi);
        }
        public Vector3d CalCFGradCorrection_bound(Vector3d[] gradPhi, double[] phi)
        {
            Vector3d aveGrad = gradPhi[OSideCell];
            return calCFCorrectPart(ref aveGrad, phi);
        }
        private double calNonOrthogonalityCorrection(Vector3d aveGrad, double[] phi)
        {
            //----------------------------修改点2-----------------------------------
#if FaceMinCor
            //!!!!!!!!以下为minCorrection
            return Vector3d.DotProduct(aveGrad, Normal - gDiff * OToN);
#else
                //!!!!!!!!!以下为overCorrection
                return Vector3d.DotProduct(aveGrad, Normal - _gT * OToN) + (phi[NSideCell] - phi[OSideCell]) * (_gT - gDiff);
#endif
        }
        public double CalNonOrthogonalityCorrection_inner(Vector3d[] gradPhi, double[] phi)//公式整合
        {
            return calNonOrthogonalityCorrection(CalAverage(gradPhi), phi);
        }
        public double CalNonOrthogonalityCorrection_bound(Vector3d[] gradPhi, double[] phi)
        {
            return calNonOrthogonalityCorrection(gradPhi[OSideCell], phi);
        }
        public double CalNonOrthogonalityCorrection_bound(out double a, out double b, Vector3d[] gradPhi, double[] phi, FvCell[] cells)//只写了匹配minCorrection的//虽然不是minCorrect也可能能用，但应该重新推导公式，进行源项线性化
        {
            a = (Math.Pow(Area, 2) - Math.Pow(gT * OToNLength, 2)) / cells[OSideCell].Volumn;
            double cor = calNonOrthogonalityCorrection(gradPhi[OSideCell], phi);
            b = cor - a * (phi[NSideCell] - phi[OSideCell]);
            return cor;
        }
        public Vector3d CalONGradient_DistanceWeighted(double[] phi)//求解Grad 梯度时使用
        {
            double deltaPhiWithDistace = (phi[NSideCell] - phi[OSideCell]) / OToNLength;
            return deltaPhiWithDistace * OToN;
        }
        public ONPair GetONCellID()
        {
            return new ONPair(OSideCell, NSideCell);
        }

    }

}
//public double calNonOrthogonalityCorrection_Old(Vector3d[] gradPhi, double[] phi) //FVM2 P244//已废弃，用新的
//{
//    return Vector3d.DotProduct(calCFCorrectedGradient_inner(gradPhi, phi), Normal)
//            - (phi[NeighborCellIndex] - phi[OwnerCellIndex]) * gDiff;
//}
